Адрес первоисточника: http://perso.wanadoo.fr/polyvalens/clemens/lifting/lifting.html

(C) C. Valens, 1999-2004

Disclaimer

This tutorial is aimed at the engineer, not the mathematician. This does not mean that there will be no mathematics, it just means that there will be no proofs in the text. In my humble opinion, mathematical papers are completely unreadable because of the proofs that clutter the text. For proofs the reader is pointed to suitable references. The equations presented are there to illustrate and to clarify things, I hope. It should not be necessary to understand all the equations in order to understand the theory. However, to understand this tutorial, a mathematical background on an engineering level is required. Also some knowledge of signal processing theory might come in handy.

The information presented in this tutorial is believed to be correct. However, no responsibilty whatsoever will be accepted for any damage whatsoever due to errors or misleading statements or whatsoever in this tutorial. Should there be anything incorrect, incomplete or not clear in this text, please let me know so that I can improve this tutorial.

This tutorial is best viewed with graphics enabled, since all special characters and equations are small .PNG files. However, care has been taken to provide alternative texts for all graphic objects.

This tutorial was developed using Microsoft's Internet Explorer 4.0, so I guess that this will be the most suitable browser to read this tutorial.

This tutorial is a sequel to the wavelet tutorial, which will fill in the blank spots on the wavelet transform map, add some detail and even explore the area outside it. We start with taking a closer look at the scaling and wavelet filters in general, what they should look like, what their constraints are and how they can be used in the inverse wavelet transform. Then we will do some algebra and develop a general framework to design filters for every possible wavelet transform. This framework was introduced by Sweldens [Swe96a] and is known as the lifting scheme or simply lifting. Using the lifting scheme we will in the end arrive at a universal discrete wavelet transform which yields only integer wavelet- and scaling coefficients instead of the usual floating point coefficients. In order to clarify the theory in this tutorial a detailed example will be presented.

In this tutorial we will go into some more detail compared to the wavelet transform tutorial, since the lifting scheme is a quite recent development and especially integer lifting [Cal96], [Uyt97b] and multi-dimensional lifting [Kov97], [Uyt97a] are not (yet) widely known. This tutorial is mainly based on [Dau97], [Cal96], [Swe96a], [Swe96b], [Cla97], [Uyt97b] and [Uyt97c].

Before we start a short note on notation. In order to be compatible with existing lifting literature this tutorial will use the same symbols. This means that the analyzing filters are denoted as and , i.e. with a tilde, while the synthesizing filters are denoted by a plain h and g. In fact, everything that has to do with the forward wavelet transform will carry a tilde. However, a slightly different form will be used here due to the fact that in [Dau97] the transposed version of is used in all calculations and because the filters and are there defined as (z^-1) and (z^-1) instead of (z) and (z). In this tutorial stands for the scaling function coefficients and for the wavelet coefficients.

As mentioned before, we can perform a wavelet transform (or subband coding or multiresolution analysis) using a filter bank. A simple one-stage filter bank is shown in figure 1 and it can be made using FIR filters. Although IIR filters can also be used, they have the disadvantage that their infinite response leads to infinite data expansion. For any practical use of an IIR filter bank the output data stream has to be cut which leads to a loss of data. In this text we will only look at FIR filters.

Figure 1
A one-stage filter bank for signal analysis and reconstruction.

figure 1 shows how a one-stage wavelet transform uses two analysis filters, a low-pass filter and a high-pass filter followed by subsampling for the forward transform. From this figure it seems only logical to construct the inverse transform by first performing an upsampling step and then to use two synthesis filters h (low-pass) and g (high-pass) to reconstruct the signal. We need filters here, because the upsampling step is done by inserting a zero in between every two samples and the filters will have to smooth this.

For the filter bank in figure 1 the conditions for perfect reconstruction are now given by [Dau97] as:

In the filter stage shown in the left part of figure 1 the signal is first filtered and then subsampled. In other words, we throw away half of the filtered samples and keep only the even-numbered samples, say. Clearly this is not efficient and therefore we would like to do the subsampling before the filtering in order to save some computing time. Let us take a closer look at what exactly is thrown away by subsampling.

Figure 2
A standard FIR filter with subsampling of the output (top) and a more efficient implementation (bottom).

In figure 2 we have drawn a standard FIR filter followed by a subsampler. If we write down its output signal y(z), just before the subsampler, for several consecutive samples, we get

.
.
.
     (2)
.
.
.

Subsampling of y(z) will now remove the middle line (row) in (1) and all the other odd lines (rows). We notice that with the odd lines removed the even-numbered filter coefficients h_e are only used with even-numbered samples x_e, while the odd-numbered filter coefficients h_o are only used with odd-numbered samples x_o. If we take the even bits together and name them h_e(z)x_e(z) and do the same with the odd bits to form h_o(z)x_o(z), then we can write the subsampled output signal y_e(z) as

     (3)

The delay z^-1 in front of the odd part in (3) comes from the delay between even and odd samples. (3) shows that we can redraw the FIR filter as shown in the right part of figure 2. In this figure we have assumed (without loss of generality) that n is even.

If we apply this remodeled FIR filter to our wavelet transform filter stage in figure 1 we end up with two equations like (3), one for each filter, so that, if we switch to vector notation, we can write

     (4)

where (z) is the polyphase matrix[1]

     (5)

The polyphase matrix now performs the wavelet transform. If we set _e(z) and _o(z) to one and we make _o(z) and _e(z) zero, i.e. (z) is the unit matrix, then the wavelet transform is referred to as the lazy wavelet transform [Swe96a]. However, I would like to rename it to the femmelet transform (femmelette being french for whimp). The femmelet transform does nothing more than splitting the input signal into even and odd components. The polyphase matrix will be used later on to build a very flexible wavelet transform.

Now let's move on to the right part of figure 1, the inverse wavelet transform. Here we have to deal with upsampling after which some filtering is performed. Upsampling is nothing more than inserting a zero in between every two samples and the consequence of this is that the filter will perform a lot of multiplications by zero, again a waste of computing time. Since the idea of moving the subsampling in front of the filters worked rather well for the forward wavelet transform, we will try a similar approach for the inverse wavelet transform, i.e. moving the upsampling behind the filters.

We look again at (2) but now imagine it as being the result of filtering an upsampled sequence of samples. If we assume that the inserted zeroes are the odd samples, then all the terms with an odd-numbered x(z) vanish and we can divide the output samples y(z) in odd and even sequences

     (6)

The "delay" z in front of the odd-numbered output samples is due to the delay between odd and even samples, necessary to merge the two sequences into the output stream y(z). As in the subsampling case we now apply this result to the reconstruction filter stage on the right in figure 1 and write down an equation similar to (4):

     (7)

     (8)

Note that when we ignore the tildes in (5), (8) is the transposed version of (5). This polyphase matrix performs the inverse wavelet transform. In case of the femmelet transform P(z) will be the unit matrix as well.

Figure 3
A one-stage filter bank for signal analysis and reconstruction using polyphase matrices.

In figure 3 we have redrawn the filter stage of figure 1, this time using the polyphase matrices. The delays we had in (4) and (7) are incorporated in the split- and merge boxes. From this figure it will be clear that the condition for perfect reconstruction now can be written as

     (9)

Here we have time-reversed one of the two polyphase matrices because in order for (9) to hold we need to cancel the delays caused by the polyphase matrices (a FIR filter is a delay line, see figure 2).

If we assume that P(z) is invertible and if we use Cramer's rule to calculate its inverse, we find

     (10)

     (11)

which implies that[2]

     (12)

     (13)

This summation is also known as a Laurent polynomial or Laurent series[3]. A Laurent polynomial differs from a normal polynomial in that it can have negative exponents. The degree of a Laurent polynomial h is defined as

     (14)

     (15)

As will be clear from our intermezzo the polyphase matrix is a matrix of Laurent polynomials and since we demanded that its determinant be equal to 1, we know that the filter pair (h,g) is complementary. The lifting theorem [Dau97] now states that any other finite filter g^new complementary to h is of the form

     (16)

where s(z²) is a Laurent polynomial. This can be seen very easily if we write g^new in polyphase form (see also [3]) and assemble the new polyphase matrix as

     (17)

As can be easily verified the determinant of the new polyphase matrix also equals 1, which proofs (16).

Similarly, we can apply the lifting theorem to create the filter ^new(z) complementary to (z) (recall the identities from (11))

     (18)

     (19)

What we just did is called primal lifting, we lifted the low-pass subband with the help of the high-pass subband. figure 4 shows the effect of primal lifting graphically.

Figure 4
Primal lifting, the low-pass subband is lifted with the help of the high-pass subband.

     (20)

     (21)

     (22)

     (23)

Dual lifting can be graphically displayed as in figure 5.

Figure 5
Dual lifting, the high-pass subband is lifted with the help of the low-pass subband.

The inverse lifting transform now also begins to take shape. If we start with a femmelet transform we only split the input stream into an even and an odd stream. Then we lift one of these streams as in the left of figures 4 and 5 by applying a Laurent polynomial to the other and adding it to the first. We can very easily undo this lifting step by again applying the same Laurent polynomial to the other stream and then subtract it from the first. In other words, inverting a lifting transform is the same as changing all the signs of the lifting Laurent polynomials in figures 4 and 5 and run it backwards, i.e. start at the output. Inverting the lifting scheme this way will always work! From this we can conclude that t(z)=-(z) and s(z)=-(z)

From the figures 4 and 5 we can see another interesting property of lifting. Every time we apply a primal or dual lifting step we add something to one stream. All the samples in the stream are replaced by new samples and at any time we need only the current streams to update sample values. In other words, the whole transform can be done in-place, without the need for auxiliary memory. This is the same as with the fast Fourier transform, where the transformed data also takes the same place as the input data. This in-place property makes the lifting wavelet transform very attractive for use in embedded applications, where memory and board space are still expensive.

We conclude this section with a note on terminology. In lifting literature the dual lifting step is also referred to as the predict step, while the primal lifting step is also referred to as the update step. The idea behind this terminology is that lifting of the high-pass subband with the low- pass subband can be seen as prediction of the odd samples from the even samples. One assumes that, especially at the first steps, consecutive samples will be highly correlated so that it should be possible to predict the odd ones from the even ones (or the other way around). The update step, i.e. lifting the low-pass subband with the high-pass subband, then is done to keep some statistical properties of the input stream, usually at least the average, of the low-pass subband.

Starting with for instance (22) we can go the other way by writing

     (24)

The form of (24) is identical to a long division with remainder of Laurent polynomials, where ^new(z) is the remainder. If we rewrite (23) a little as well, we obtain

     (25)

and we see that for the polyphase matrix we have to perform two long divisions in order to extract one lifting step. Once we have extracted one such lifting step we can continue by extracting more lifting steps from the new polyphase matrix until we end up with the unit matrix, or a matrix with only two constants on its main diagonal. In fact, in [Dau97] it is proven that it will always be possible to do this when we start with a complementary filter pair (h,g), i.e. P(z) can always be factored into lifting steps:

     (26)

In (26) K₁ and K₂ are scaling constants unequal to zero. If scaling is not desired for some reason, it is even possible to factor the scaling matrix into four more lifting steps, one of which can be combined with the last real lifting step so that factoring a scaling matrix costs three extra lifting steps [Dau97].

To perform a long division with remainder on Laurent polynomials we can use the Euclidean algorithm for Laurent polynomials[4]. The Euclidean algorithm was originally developed to find the greatest common divisor (gcd) of two natural numbers, but we can extend it to Laurent polynomials as well. In [Dau97] it is given as:

Take two Laurent polynomials a(z) and b(z) 0 with |a(z)| |b(z)|. Let a₀(z) = a(z) and b₀(z) = b(z) and iterate the following steps, starting from i = 0

     (27)

The `%'-symbol in the second line of (27) means mod, i.e. integer divide with remainder but only keeping the remainder, and it is the same symbol as used in the C programming language for the mod operation. The `/'-symbol in the third line is the C-language div operator. The result of this algorithm can be written as

     (28)

     (29)

     (30)

     (31)

     (32)

     (33)

With the help of (11) we can now assemble the synthesizing polyphase matrix as well. First we check the determinant of (33):

     (34)

so we have to scale (11) a bit before using the equalities to get:

     (35)

Do not get confused here, P(z) is not time-reversed! Remember, if we want to check (9) we have to time-reverse one of the polyphase matrices. From (35) we can find the synthesizing filters as follows:

     (36)

     (37)

     (38)

The next step is the factoring of the polyphase matrices into lifting steps. We start with the extraction of a dual lifting step[5].

     (39)

     (40)

     (41)

Note also that we have found three greatest common divisors. If we choose the middle line of (41) as the factorization we have a symmetrical one, which goes nicely with (z) as well, and we arrive at the following decomposition:

     (42)

We can continue by extracting a primal lifting step. For this we apply the Euclidean algorithm to _e(z) and _o(z) of (42), almost not worth mentioning it, and find:

     (43)

This equation gives a fully factored version of the filters from (29). If we finally use (43) with (4) we can display our wavelet transform graphically as in figure 6 while the corresponding wavelet and scaling function are displayed in figure 7.

Figure 6
The implementation of the wavelet transform of this example.

From figure 6 we can generalize the lifting steps as:

     (44)

Figure 7
The scaling function (left) and the wavelet (right) that go with this example.

The lifting scheme has some properties which are not found in many other transforms. figure 6 shows a few of these properties and we will now discuss some of the most important.

The inverse transform is immediately clear: change the signs of all the scaling factors, replace "split" by "merge" and go from right to left, i.e. reverse the data flow. This easy invertibility is always true for the lifting scheme.

Lifting can be done in-place (44): we never need samples other than the output of the previous lifting step and therefore we can replace the old stream by the new stream at every summation point. Not immediately clear from this figure is that when we iterate a filter bank using in-place lifted filters we end up with interlaced coefficients. This can be seen as follows. We split the input in odd- and even-numbered samples and perform the in-place lifting steps. After one complete step the high-pass filtered samples, the wavelet coefficients, sit in the odd-numbered places and the low-pass filtered samples sit in the even-numbered places. Then we perform another transform step, but only using the low-pass filtered samples, so that this sequence will again be divided into odd- and even-numbered samples. Again the odd-numbered samples are transformed into wavelet coefficients, while the even-numbered samples will be processed further so that in the end all wavelet coefficients will be interlaced.

The third important property has not been mentioned yet, but it shows clearly from figure 6: lifting is not causal. Usually this is not really a problem, we can always delay the signal enough to make it causal, but it will never be real- time. In some cases however it is possible to design a causal lifting transform.

The last important property we will mention here is the calculation complexity. In [Dau97] it is proven that for long filters the lifting scheme cuts computation complexity in half, compared to the standard iterated FIR filter bank algorithm. This type of wavelet transform has already a complexity of N, in other words, much more efficient than the FFT with its complexity of Nlog(N) and lifting speeds things up with another factor of two. This is where the title of this tutorial comes from: it is a fast wavelet transform and therefore we will refer to it as the fast lifting wavelet transform of FLWT.

The last stage of our voyage to the ultimate[6] wavelet transform is the stage where we make sure that the wavelet coefficients are integers. In classical transforms, including the non-lifted wavelet transforms, the wavelet coefficients are assumed to be floating point numbers. This is due to the filter coefficients used in the transform filters, which are usually floating point numbers. In the lifting scheme it is however rather easy to maintain integer data, although the dynamic range of the data might increase. That this is possible in the lifting scheme has to do with the easy invertibility property of lifting.

The basic lifting step is given in (44) and we rewrite it here a little modified as [Uyt97c]:

     (45)

     (46)

where we have used to denote the rounding operation. (46) is fully reversible:

     (47)

We have to take care however, because we did not consider the scaling step in the previous paragraph. Scaling usually does not yield integer results but it is a part of the lifting transform. The simplest solution to this problem is to forget all about scaling and just keep in mind that the transform coefficients actually have to be scaled. This is important for instance in denoising applications. If scaling is ignored, then it is desirable to let the scaling factor be as close to one as possible. This can be done using the non-uniqueness of the lifting factorization. Another solution is to factor the scaling into lifting steps as well [Dau97].

As mentioned before the integer lifting transform can not guarantee the preservation of the dynamic range of the input signal. Usually the dynamic range doubles [Uyt97c], but there are schemes that can keep the dynamic range. In [Cha96] an interesting lifting transform with the so-called property of precision preservation (PPP) is described. This transform makes use of the two-complement representation of integers in a computer and the wrap- around overflows cause in this representation. The disadvantage of such a transform is that large coefficients may be represented by small values and it is therefore difficult to take decisions on coefficient values.

We have now finished our self-imposed task of transforming the CWT into a practical implementation. In this tutorial we have seen how we can use the lifting scheme to build a very versatile wavelet transform. After first optimizing the subsampled and upsampled FIR filters from the wavelet tutorial, through the use of some algebra we arrived at a scheme to build a wavelet transform using primal and dual lifting blocks. These modules allowed us to build any wavelet transform, which fits in the classical framework, and more. Adapting the lifting scheme we will be well armed: amongst our weaponry are such elements as[7] easy invertibility of any transform, in-place calculation of the transform and easy integer transform coefficients without losing any of its features. And there are many more features [Cal96], [Uyt97b], [Dau97].

However, this does not mean that this is the only way to go. There are probably as many wavelet transforms as there are wavelets. Due to the infinite variety of wavelets it is possible to design a transform which maximally exploits the properties of a specific wavelet[8], and of course this has been done. While researching wavelet theory I have stumbled upon morlets[9], coiflets, wavelants, slantlets, brushlets and wavelet packets to name a few. The lifting scheme on the other hand is a really general scheme, which makes it very suitable for experimenting while the in-place and integer properties make it extremely useful for embedded systems where memory is still expensive. With the application described in this report in mind, it will be clear that these are the reasons for studying the lifting scheme.

Finally, four remarks to conclude this tutorial:

• In [Uyt97b] the filter factoring algorithm is used to split the original filters in simpler filters and one primal and dual lifting step. These lifting steps are then used to make the original wavelet transform integer. This is some kind of hybrid (trans)form but very effective.
• Up to now we have only spoken about one-dimensional transforms. It is however easily possible to extend the lifting transform to the multi-dimensional case. Not only can the lifting transform be used in a classical separable multi- dimensional setting, but it can be made truly multi-dimensional. In [Uyt97a] the lifting transform is extended to a true two-dimensional transform, while in [Kov97] the complete theoretical foundations are laid out for any dimension. The principles behind lifting do not change at all in the multi-dimensional setting.
• One of the advantages of the lifting scheme as pointed out in for instance [Dau97] and [Cal96] is that the lifting scheme allows for an introduction into wavelet theory without the use of Fourier theory. We do not agree with this on the grounds that from the lifting scheme it is totally unclear why there should be wavelets in it at all. The concept of wavelets is completely unnecessary to understand the lifting scheme and therefore, we feel, it should not be used as an introduction to wavelet theory.
• The lifting scheme is constantly under development and is investigated by many. Recent additions are the lifting scheme in a redundant setting in order to improve the translation invariance [Sto98] and adaptive prediction schemes for integer lifting [Cla97].

11. Notes

12. References